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1.  Introduction 


The  behavior  of  materials  under  conditions  of  extreme  temperature  and  pressure  is  of  critical 
importance  in  many  fields  of  physics  and  fluid  science  (1-3).  The  properties  of  materials  at 
these  conditions  can  be  measured  through  shock  experiments,  which  are  capable  of  producing 
pressures  up  to  several  hundred  GPa  and  temperatures  exceeding  10,000  K.  Moreover, 
determining  the  shock  properties  of  energetic  materials  is  a  crucial  task  in  the  field  of  detonation 
science  (4). 

Experimental  measurements  of  the  properties  of  shocked  materials  are  often  difficult  because 
instrumentation  must  be  capable  of  spanning  a  wide  range  of  pressures  (1-300  GPa)  and 
temperatures  (500-15000  K).  Additionally,  energy  releases  that  might  accompany  a  material 
under  shock  conditions,  as  well  as  the  time  and  length  scales  over  which  the  event  occurs,  have 
thwarted  extensive  experimental  studies  of  many  fundamental  substances.  As  a  further 
complication,  recent  theoretical  predictions  suggest  that  the  detonation  products  of  some  systems 
supercritical  phase  separate,  significantly  altering  the  shock  properties  (e.g.,  references  [5-7]. 
Unfortunately,  current  experimental  techniques  are  not  capable  of  delineating  the  phase 
separation  of  materials  under  shock  conditions;  thus,  this  behavior  has  yet  to  be  verified. 

Such  laboratory  challenges  have  necessitated  the  development  of  theoretical  predictive 
capabilities  to  complement  the  experimental  analyses.  To  date,  the  most  reliable  theoretical 
treatments  for  predicting  shock  properties  apply  statistical  mechanical  approaches  such  as 
variational  perturbation  theory  (e.g.,  Ross  [5]  or  integral  equation  theory  (e.g.,  references 
[9-11]).  These  approaches  predict  the  shock  properties  by  minimizing  the  Gibbs  free-energy  and 
by  requiring  that  the  total  number  of  elements  constituting  the  chemically  reacting  species  is 
conserved.  Thermochemical  software  such  as  the  chemical  equilibrium  code  (12)  and  Cheetah 
(13)  are  capable  of  performing  such  calculations.  However,  these  approaches  require  accurate 
equations  of  state  (EOS)  for  the  reactive  mixture,  thus  limiting  the  prediction  of  the  shock 
properties  of  materials  to  those  for  which  equations  of  state  are  available.  The  calculations  also 
require  heats  of  formation  and  densities  of  the  unshocked  material.  These  quantities  are  often 
unknown  for  notional  or  novel  energetic  materials  as  well  as  for  materials  in  highly  nonideal 
environments  (e.g.,  energetic  materials  packed  in  polymer  matrices  or  confined  in  carbon 
nanotubes).  Furthermore,  approximations  must  often  be  made  within  the  theoretical  models  to 
keep  the  calculations  tractable,  for  example,  to  develop  an  analytical  representation  of  the  fluid 
equation  of  state  (11)  or  in  applying  the  van  der  Waals  one-fluid  approximation  (12, 13).  These 
types  of  approximations  can  add  uncertainty  to  the  predictive  capabilities  of  the  methods. 

A  powerful  simulation  method  available  for  studying  the  shock  properties  of  materials  while 
providing  insight  into  atomic-level  phenomena  is  the  molecular  dynamics  (MD)  method  (14-31). 
The  method  can  be  used  irrespective  of  rate  limitations,  the  production  of  huge  energy  releases, 
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extreme  thermodynamic  conditions,  or  other  regimes  that  are  presently  inaccessible  by 
experimental  methods.  MD  evaluation  of  the  Hugoniot  states  of  a  material  can  be  accomplished 
by  calculating  properties  behind  the  shock  discontinuity  in  a  shockwave  simulation  (30)  or  by 
generating  an  EOS  for  subsequent  evaluation  of  the  Hugoniot  conservation  relations  (24,  27). 
Recently,  an  equilibrium  MD  method  has  been  introduced,  termed  uniaxial  Hugoniostat  (29), 
which  utilizes  equations  of  motion  that  constrain  the  system  during  the  MD  simulation  such  that 
the  time-averaged  properties  correspond  to  those  on  the  shock  Hugoniot  curve.  For  evaluating 
the  shock  Hugoniot  of  a  material,  the  uniaxial  Hugoniostat  technique  is  more  efficient  than  the 
method  of  generating  an  EOS  using  standard  MD  for  subsequent  evaluation  of  the  conservation 
equations  (24,  27).  A  significant  drawback  of  all  MD  approaches,  however,  is  that  they  require 
an  accurate  model  of  the  interaction  potential  experienced  between  all  species  in  the  shocked  and 
unshocked  states.  If  the  relative  species  concentrations  of  the  products  in  the  shocked  state  of 
interest  are  not  known,  then  the  MD  approach  requires  an  interaction  potential  that  simulates 
bond  breaking  and  bond  formation  in  order  to  establish  the  chemical  equilibrium  of  the  final 
shocked  state.  Although  significant  advances  have  been  made  in  developing  potentials  that 
reproduce  the  characteristics  of  a  detonation,  the  potentials  are  still  highly  idealized 
representations  of  the  energetic  molecular  system  (14-23,  25-28).  However,  if  the  relative 
species  concentrations  are  known  at  the  conditions  of  the  shocked  state,  then  generating  the 
shock  Hugoniot  curve  using  either  the  conventional  or  the  uniaxial  Hugoniostat  MD  method  only 
requires  appropriate  potentials  that  describe  nonreacting  interactions  among  all  species  present. 

One  of  several  alternative  methods  for  calculating  chemically  reactive  systems  which  is 
appropriate  for  determining  the  shock  properties  of  materials  is  Reaction  Ensemble  Monte  Carlo 
(RxMC)  (32,  33).  RxMC  circumvents  some  of  the  problems  associated  with  conventional  and 
uniaxial  Hugoniostat  MD  methods.  RxMC  requires  neither  a  priori  knowledge  of  the  relative 
concentrations  of  the  species  in  the  shocked  state  nor  a  potential  that  describes  bond  breaking  or 
bond  formation.  The  method  only  requires  (a)  functions  that  accurately  describe  nonreactive 
interactions  between  all  possible  species  in  the  equilibrium  mixture  of  the  final  shocked  state 
(i.e.,  intermolecular  potentials)  and  (b)  ideal  gas  partition  functions  for  all  species  in  the  mixture. 
For  a  given  intermolecular  potential  model,  the  method  will  provide  infonnation  on  the  chemical 
equilibrium  state,  such  as  the  density  of  the  reactive  mixture,  the  mole  fractions  of  reactive 
species,  the  change  in  the  total  number  of  moles,  and  the  internal  energy.  The  intermolecular 
model  can  contain  various  levels  of  detail  including  multisite  molecules  and  electrostatic 
contributions.  Numerous  reactions  can  be  simulated  simultaneously  in  multiple  phase  systems 
by  perfonning  Monte  Carlo  sampling  of  forward  and  reverse  reaction  steps.  Using  the  RxMC 
method,  the  dependence  of  the  chemical  equilibria  on  system  conditions  such  as  temperature, 
pressure,  and  the  surrounding  environment  (e.g.,  a  condensed  phase  or  highly  nonideal 
environment)  can  be  studied. 
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In  this  work,  we  demonstrate  the  applicability  of  RxMC  for  calculating  the  shock  Hugoniot 
properties  of  materials.  We  illustrate  the  method  on  shocked  liquid  N2,  which  dissociates  into 
atomic  nitrogen  (N)  within  a  particular  thermodynamic  regime.  This  reaction  has  been  studied 
extensively  by  both  experimental  and  theoretical  techniques.  This  report  is  the  first  of  two 
reports  that  illustrate  the  RxMC  method.  In  the  second  report  (34),  we  consider  shocked  liquid 
NO,  which  is  (nearly)  an  irreversible  decomposition  reaction  that  generates  a  mixture  of 
homonuclear  products  (N2  and  02). 

The  outline  of  the  report  is  as  follows.  The  RxMC  methodology  applied  to  the  simulation  of  the 
shock  properties  of  materials  is  described  in  section  2.  Simulation  details  and  models  can  be 
found  in  section  3  and  application  of  the  method  to  shocked  liquid  N2  is  presented  in  section  4. 


2.  Methodology 


2.1  Reaction  Ensemble  Monte  Carlo 

The  RxMC  method  (32,  33)  is  designed  to  minimize  the  Gibbs  free-energy,  thus  detennining  the 
true  chemical  equilibrium  state  irrespective  of  rate  limitations.  RxMC  requires  intennolecular 
potentials  for  the  molecular  species  that  are  present  in  the  reactive  mixture  and  often  uses 
spherically-averaged  potentials  such  as  Lennard- Jones  or  Exponential-6  models  (35).  RxMC 
also  requires  inputting  the  ideal-gas  internal  modes  (vibration,  rotation,  and  electronic)  for  each 
reactive  species.  These  contributions  can  be  included  by  calculating  internal  partition  functions 
from  molecular  energy-level  data  (32)  or  by  using  tabulated  thermochemical  data  (33). 
Regardless  of  the  approach  taken,  the  required  information  is  readily  available  in  standard 
sources  (35-37)  or  can  be  generated  using  quantum  mechanical  calculations.  Finally,  the 
particular  reactions  occurring  in  the  system  must  be  specified.  Provided  that  a  sufficient  set  of 
independent  reactions  are  specified,  this  requirement  is  not  a  considerably  limiting  factor  since 
insignificant  reactions  are  easily  discernable  by  negligible  product  concentrations. 

RxMC  implementation  provides  information  on  the  chemical  equilibrium  state,  such  as  the 
density  of  the  reactive  mixture,  mole  fractions  of  reactive  species,  the  change  in  the  total  number 
of  moles,  and  the  internal  energy.  RxMC  can  be  performed  in  many  different  types  of 
ensembles,  including  canonical,  isothermal-isobaric,  Gibbs  (38),  and  other  less  common 
ensembles  (39).  Furthennore,  RxMC  can  be  performed  for  multiple  reactions  and  multiple 
phases  (32,  33,  38,  40-46).  RxMC  does  not  simulate  bond  breaking  or  forming;  these  relatively 
rare  events  in  standard  Metropolis  Monte  Carlo  (47)  would  result  in  considerable  statistical 
uncertainty.  Rather,  RxMC  directly  samples  forward  and  reverse  reaction  steps  as  Monte  Carlo- 
type  moves  according  to  the  stoichiometry  of  the  reactions  being  sampled.  The  isothermal- 
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isobaric  version  of  the  RxMC  method  containing  J  number  of  reactions  involves  the  following 
trial  moves: 

1 .  a  change  in  the  position  or  orientation  of  a  molecule,  chosen  at  random; 

2.  a  forward  step  for  randomly  chosen  reaction  j,  in  which  reactant  molecules  are  chosen  at 
random  and  changed  to  product  molecules; 

3.  a  reverse  step  for  randomly  chosen  reaction  j,  in  which  product  molecules  are  chosen  at 
random  and  changed  to  reactant  molecules;  and 

4.  a  random  change  in  the  simulation  box  volume. 

Step  1  ensures  that  thennal  equilibrium  is  established  for  the  user-specified  temperature, 
steps  2  and  3  ensure  that  chemical  equilibrium  is  established,  while  step  4  satisfies  the 
requirement  of  mechanical  equilibrium  for  the  user-specified  pressure.  Note  that  the  transition 
probabilities  in  steps  2  and  3  do  not  require  specifying  the  values  of  the  chemical  potentials  nor 
chemical  potential  differences  for  any  of  the  mixture  components. 

The  acceptance  probability  of  state  k  going  to  state  /  for  a  particle  displacement  or  orientation 
step  is  given  by 

P?  =min{l,exp(-/?A<y„)},  (1) 


where  AC/*/  =  U/-  Uk  is  the  change  in  the  configurational  energy  and  ft=  1/k B7)  kB  is  the 
Boltzmann  constant  and  T  is  the  temperature. 


The  acceptance  probability  for  a  reaction  step  is  given  by 


P%  =  min^ 


\ 

> 

v]itj 

avn  ( /?/\  T  I  ) 

1  A3,  J 

expt  u  ki ) 

(2) 


where  Cj  is  the  total  number  of  species  in  reaction  j;  Nj  is  the  total  number  of  molecules  of  species 
i;  Vp  is  the  stoichiometric  coefficient  of  species  i  in  reaction  /;  c,  is  the  molecular  extent  of 
reaction  for  reaction  j;  qmui  is  the  quantum  partition  function  for  the  internal  modes  of  an  isolated 
molecule  of  species  i,  which  includes  vibrational,  rotational,  and  electronic;  A,  is  the  thermal  de 
Broglie  wavelength  of  species  i;  and  V  is  the  total  volume  of  the  system.  Equation  2  is 
appropriate  for  both  forward  and  reverse  reaction  steps  (£,/  =  1  for  a  forward  step  and  cy  =  - 1  for  a 
reverse  step),  where  the  stoichiometric  coefficients  are  taken  to  be  positive  for  product  species 
and  negative  for  reactant  species.  (For  example,  consider  the  reaction  A<=>2B.  For  the  reaction 
as  written,  the  stoichiometric  coefficients  for  species  A  and  B  are  vA=  -1  and  vB  =  +2  while 
£,  =  +1.  Then  for  the  reverse  reaction  step,  vA  and  vB  again  are  -1  and  +2,  respectively,  but  now 

S=-l.) 
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Finally,  a  random  change  in  the  simulation  box  volume  is  accepted  with  the  probability 


Plf  =  mini  1,  exp 


„(y,-y,)+N  in  T 


V, 


kj 


(3) 


where  Pimp  is  the  user-specified,  or  imposed,  pressure.  Derivations  of  these  transition 
probabilities  along  with  further  details  of  the  methodology  can  be  found  in  the  original  papers 
(32,  33,  38). 

2.2  Calculation  of  Shock  Hugoniot  Properties 

The  thermodynamic  quantities  of  a  material  in  the  initial  unshocked  state  and  the  final  shocked 
state  are  related  by  the  conservation  equations  of  mass,  momentum,  and  energy  across  the  shock 
front  as  follows  ( 4 ): 


Mass: 

ll 

B 

i 

(4) 

Momentum:  P  —  P0  =  p0uD. 

(5) 

Energy: 

E-E0  =  y2(P  +  p0xv0-V) . 

(6) 

In  equations  4-6,  E  is  the  specific  internal  energy,  P  is  the  pressure,  p  is  the  specific  density, 

V  =  Up  is  the  specific  volume,  D  is  the  velocity  of  the  shock  wave  propagating  through  the 
material,  and  u  is  the  mass  velocity  of  the  products  behind  the  shock  wave.  The  term  “specific” 
refers  to  the  quantity  per  unit  mass,  while  the  subscript  “o”  refers  to  the  quantity  in  the  initial 
unshocked  state. 

The  shock  wave  velocity  D  can  be  calculated  by  solving  equations  4  and  5  for  the  mass  velocity 
u  and  equating  these  expressions.  The  resulting  expression,  termed  the  Rayleigh  line,  can  be 
written  as 

R  =  p„2D 2  -  (P  -  P0)(Vo  -  V)  =  0  .  (7) 

The  so-called  Hugoniot  function  satisfies  equation  6  as 

Hg (T,V)  =  0  =  E-E0-  Vi{P  +  Po)(V0  -  V) .  (8) 

Note  that  the  quantities  E  and  V  are  extensive  quantities  in  equations  6-8  and  thus  dependent  on 
the  relative  amounts  of  the  reactive  species.  The  extensive  quantities  used  in  this  work  were 
formulated  on  a  specific  basis  (per  gram);  alternatively,  these  quantities  can  be  fonnulated  on  a 
total  system  basis  (total  number  of  moles).  The  relative  amounts  of  the  reactive  species  along 
with  the  quantities  E,  P,  and  V  are  calculated  explicitly  during  the  RxMC  simulation.  The 
internal  energy  is  calculated  during  the  RxMC  simulation  based  upon  the  following  derivation. 
The  thermodynamic  definition  of  the  internal  energy  is  given  as 

E  =  H -PV ,  (9) 
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where  H  is  the  specific  enthalpy  that  can  be  written  as  a  sum  of  the  ideal-gas  (H°)  and  excess 
enthalpies  (//c): 

H  =  H°+He.  (10) 

The  ideal-gas  and  excess  enthalpies  can  be  written  as 

»■=&,»;  an 

1=1 


and 

He  =Ucont  +  PV-RT ,  (12) 

so  that 

cj 

H  =  YjyiH°  +Uconf  +PV-RT.  (13) 

i= 1 

lfont  is  the  total  configurational  energy  calculated  during  the  simulation  from  the  species-species 
interactions;  y,  is  the  mole  fraction  of  species  i;  Cj  is  the  total  number  of  species;  H°  is  the 
specific  ideal-gas  enthalpy  of  pure  species  i,  which  can  be  determined  solely  from  tabulated 
thermochemical  data  at  the  appropriate  temperature,  T  (37,  48,  49),  or  with  tabulated 
thennochemical  data  supplemented  with  computed  values  where  data  is  lacking  (e.g.,  [50,  5/]); 
and  R  is  the  universal  gas  constant.  Substituting  equation  13  into  equation  9, 

ci 

E  =  YjyiH^+  jjconf  -R  T  ,  (14) 

i= 1 

which  is  the  specific  internal  energy  expression  needed  in  equation  8. 

The  search  algorithm  for  locating  a  point  on  the  Hugoniot  curve  used  in  this  study  is  as  follows: 

Step  1 :  For  a  user-specified  pressure,  P,  a  few  RxMC-NPT  simulations  were  performed  at 
temperatures  believed  to  be  near  Hg(T,V)  =  0. 

Step  2:  A  functional  form  for  the  Hg  vs.  T  plot  is  determined,  e.g.,  fitted  to  a  quadratic 
polynomial,  and  the  temperature  (THg)  that  satisfies  equation  8  at  Hg  =  0  is 
interpolated. 

Step  3:  Similarly,  functional  forms  for  the  other  desired  quantities  (V,  E,  D)  were  determined 
and  interpolated  at  THg. 

Depending  on  the  initial  guess  of  THg,  additional  RxMC  simulations  may  be  required  to  achieve 
the  desired  accuracy.  Typically,  4-6  RxMC  simulations  are  needed  to  detennine  a  single  point 
on  the  Hugoniot  curve,  where  best  results  are  obtained  when  chosen  values  of  T  include  both 
Hg>0  and  Hg<0.  Steps  1-3  are  then  repeated  to  trace  out  the  entire  Hugoniot  curve. 
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3.  Simulation  Model  and  Details 


3.1  Intermolecular  Potential  Models 

The  species  particles  interact  through  the  exponential-six  potential,  which  can  be  expressed  as: 


U  , (r)  = 

exp -6  v  ^ 


e 

6 

—  exp 

( 

a 

1-^- 

> 

ro 

6 

1 

_ 1 

a 

V 

Tm 

J 

V  r  J 

- 

r  <  r 


r  >  r 


(15) 


where  s  is  the  depth  of  the  attractive  well  between  particles,  rm  is  the  radial  distance  at  which  the 
potential  is  a  minimum,  while  a  controls  the  steepness  of  the  repulsive  interaction.  The  cutoff 
distance  rcore  is  included  to  avoid  the  unphysical  singularity  in  the  potential  function  as  r— >0. 

The  potential  parameters  for  the  species  considered  in  this  work  are  given  in  table  1 .  A  spherical 
cutoff  of  2.5rmjN2  was  applied  with  long  range  corrections  added  to  account  for  interactions 
beyond  this  distance  (52).  Electrostatic  contributions  were  ignored  between  species.  The  unlike 
interactions  between  species  i  and  j  were  approximated  by  the  Lorentz-Berthelot  mixing  rules 
(55)  for  Sjj,  (Xij,  and  rm4j\ 

Sij  (£i£j)  •>  Gtjj  (oCiOCj)  ”  ,  r mjj  (fmj  +  rmj)/2  , 


while 


r  core, zy 


=  (r, 


core,;  T  core^')/2  . 


Table  1.  Exponential-6  potential  parameters. 


(16) 


Species 

f  core 

(A) 

f  lit 

(A) 

e/kB 

(K) 

a 

Source 

(Reference  No.) 

n2 

1.13 

4.2005 

101.10 

12.684 

(ID 

N 

0.98 

2.5688 

88.181 

11.013 

(ID 

The  vibrational  and  rotational  contributions  to  the  ideal-gas  partition  function  of  N2  used  in 
simulating  the  N2  dissociation  reaction  were  calculated  using  a  standard  source  (36)  and 
supplemented  with  electronic  level  constants  that  included  the  ground  state  and  six  excited 
electronic  states  (48).  For  N,  the  electronic  energy  levels  were  taken  from  Moore  and  Gallagher 
(49).  The  corresponding  thermochemical  reference  data  were  used  in  calculating  the  ideal-gas 
enthalpies  (H°)  required  in  equation  11  (57,  48,  49). 
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3.2  Simulation  Details 


Constant-pressure  RxMC  simulations  of  shocked  N2  were  initiated  from  3375  N2  particles  placed 
on  a  face-centered-cubic  lattice  structure.  The  standard  periodic  boundary  conditions  and 
minimum  image  convention  were  used  (53).  Simulations  were  performed  in  steps,  where  a  step 
(chosen  with  equal  probability)  was  either  a  particle  displacement,  forward  reaction  step,  or 
reverse  reaction  step.  A  change  in  the  simulation  cell  volume  was  attempted  every  2500  steps. 
Simulations  were  equilibrated  for  0.3  x  107  steps  after  which  averages  of  the  quantities  were 
taken  over  2.0  x  107  steps.  Uncertainties  were  estimated  using  the  method  of  block  averages  by 
dividing  the  production  run  into  10  equal  blocks  (52).  Reported  uncertainties  are  one  standard 
deviation  of  the  block  averages.  The  maximum  displacement  and  volume  change  were  adjusted 
to  achieve  an  acceptance  fraction  of  -0.33  and  0.5,  respectively.  Depending  on  the  system 
conditions,  the  acceptance  fraction  of  the  reaction  steps  ranged  from  0.075-0.375.  Calculated 
quantities  were  reduced  by  the  exponential-6  potential  energy  ( s )  and  size  (rm)  parameters  of  N2. 


4.  Application 


We  consider  shocked  N2  in  the  pressure  range  3-90  GPa  for  which  reliable  experimental  data  are 
available  (54,  55).  At  pressures  between  -30-100  GPa,  the  N=N  triple  bond  is  destabilized  and 
molecular  nitrogen  dissociates  into  atomic  nitrogen  (56,  57).  At  higher  pressures,  theory 
suggests  that  nitrogen  can  exist  as  a  metastable  polymeric  phase  of  N  atom  clusters  that  are 
covalently  bonded  (58-62),  before  losing  this  covalency  at  still  higher  pressures.  In  the  present 
work,  we  consider  only  the  regime  where  molecular  nitrogen  is  believed  to  dissociate  into  atomic 
nitrogen,  N2<=>2N.  Analogous  to  the  work  of  Fried  and  Howard  (11),  we  consider  two  models 
for  this  reaction,  a  reactive  model  that  includes  the  dissociation  reaction,  and  a  nonreactive 
model  that  does  not. 

We  determined  the  shock  Hugoniot  properties  of  liquid  N2  using  the  calculated  initial  states 
given  in  table  2.  The  values  given  in  table  2  were  detennined  by  perfonning  a  canonical 
ensemble  (constant-N  VT)  Monte  Carlo  simulation  of  N  =  3375  N2  molecules  at  T  =  77.0  K  and 
at  a  specific  volume  of  V  =  1.238  cm  /g.  A  comparison  of  the  pressure  and  internal  energy 
determined  from  this  NUT  simulation  with  experiment  is  given  in  table  2.  The  shock  Hugoniot 
properties  were  determined  by  carrying  out  the  prescription  outlined  in  section  2.2.  The  raw  data 
determined  from  a  series  of  constant-pressure  RxMC  simulations  at  several  different 
temperatures  are  given  in  table  3.  Also  reported  in  table  3  are  the  values  of  the  Hugoniot 
expression  (equation  8)  using  the  predicted  thennodynamic  data.  Quadratic  polynomials  were 
used  in  the  fitting  procedure  of  steps  2  and  3  (see  section  2.2)  with  the  exception  of  the  shock 
wave  velocity  (D),  where  a  linear  equation  was  used.  A  comparison  of  the  shock  properties 
along  the  principal  Hugoniot  calculated  from  the  RxMC  simulations  and  the  available 
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Table  2.  Initial  fluid  state  used  to  evaluate  equation  8. 


Thermodynamic 

Property 

Liquid  N2 

Experiment  (55) 

N 1T-MC 

Temperature,  T  (K) 

77.0  ±0.5 

77.0 

Density,  p  (g/cm3) 

0.808  ±0.003 

0.808 

Pressure,  P  (MPa) 

— 

50.49  ±0.02 

Energy,  E  (kJ/g) 

— 

-0.441  ±0.004 

Table  3.  Constant-pressure  RxMC  simulations  of  shocked  liquid  N2. 


T 

<P> 

Mole  Fraction0 

<V> 

<C7conf> 

H°b 

Ec 

Hgd 

Dc 

(K) 

(GPa) 

<x(N2)> 

<  x(N)> 

(cm3/g) 

(kJ/g) 

(kJ/g) 

(kJ/g) 

(kJ/g) 

(km/s) 

Pimp  =  2.96  GPa 

475 

3.010(2) 

1.0000 

0.0000 

0.7532(8) 

0.182(2) 

0.1852 

0.227 

-0.0732 

3.059 

500 

3.008(1) 

1.0000 

0.0000 

0.7572(7) 

0.188(2) 

0.2110 

0.250 

-0.0429 

3.071 

525 

3.008(2) 

1.0000 

0.0000 

0.7606(9) 

0.192(3) 

0.2379 

0.274 

-0.0138 

3.082 

550 

3.008(2) 

1.0000 

0.0000 

0.7646(9) 

0.197(3) 

0.2643 

0.298 

0.0165 

3.094 

Pimp  =  4.74  GPa 

850 

4.795(3) 

1.0000 

0.0000 

0.7047(6) 

0.463(4) 

0.5955 

0.806 

-0.0433 

3.693 

900 

4.795(3) 

0.9998(1) 

0.0002(1) 

0.7090(8) 

0.471(3) 

0.6566 

0.860 

0.0210 

3.708 

950 

4.790(3) 

0.9997(1) 

0.0003(1) 

0.7132(6) 

0.478(5) 

0.7201 

0.916 

0.0879 

3.721 

Pimp  =  10.0  GPa 

1950 

10.072(7) 

1.0000 

0.0000 

0.6197(9) 

1.249(7) 

1.9407 

2.611 

-0.0751 

4.984 

2000 

10.081(4) 

0.9998(1) 

0.0002(1) 

0.6210(6) 

1.254(9) 

2.0101 

2.670 

-0.0115 

4.992 

2050 

10.066(8) 

0.9998(1) 

0.0002(1) 

0.6235(6) 

1.260(9) 

2.0749 

2.726 

0.0609 

4.998 

Pimp  =  18.1  GPa 

3850 

18.195(2) 

0.9997(1) 

0.0003(1) 

0.5549(7) 

2.398(2) 

4.4512 

5.707 

-0.0806 

6.380 

3900 

18.191(2) 

0.9997(1) 

0.0003(1) 

0.5557(7) 

2.401(1) 

4.5179 

5.761 

-0.0168 

6.383 

3950 

18.189(2) 

0.9997(1) 

0.0003(1) 

0.5569(8) 

2.407(2) 

4.5851 

5.820 

0.0533 

6.388 

Pimp  =  29.9  GPa 

6700 

30.023(2) 

0.9930(1) 

0.0070(1) 

0.5001(8) 

3.963(1) 

8.5267 

10.502 

-0.1103 

7.890 

6725 

30.015(2) 

0.9939(1) 

0.0061(1) 

0.5017(8) 

3.960(2) 

8.5392 

10.503 

-0.0810 

7.897 

6750 

30.008(3) 

0.9936(1) 

0.0064(1) 

0.5024(7) 

3.966(3) 

8.5808 

10.543 

-0.0274 

7.900 

6775 

30.009(2) 

0.9926(1) 

0.0075(1) 

0.5008(7) 

3.954(2) 

8.6431 

10.586 

-0.0089 

7.892 

6800 

30.022(2) 

0.9923(1) 

0.0077(1) 

0.5010(7) 

3.957(1) 

8.6847 

10.624 

0.0271 

7.895 

Pimp  =  36.0  GPa 

7900 

36.083(3) 

0.9762(1) 

0.0238(1) 

0.4762(9) 

4.666(2) 

10.8830 

13.204 

-0.1105 

8.514 

7950 

36.121(3) 

0.9753(1) 

0.0247(1) 

0.4765(6) 

4.672(2) 

10.9868 

13.299 

-0.0252 

8.520 

7975 

36.100(2) 

0.9750(1) 

0.0250(1) 

0.4766(5) 

4.668(2) 

11.0363 

13.338 

0.0233 

8.518 

8000 

36.143(2) 

0.9744(1) 

0.0256(1) 

0.4766(5) 

4.677(2) 

11.0924 

13.395 

0.0644 

8.523 

8050 

36.141(2) 

0.9734(1) 

0.0266(1) 

0.4768(7) 

4.675(2) 

11.2008 

13.486 

0.1606 

8.524 

8475 

36.120(3) 

0.9648(1) 

0.0352(1) 

0.4796(11) 

4.686(2) 

12.1361 

14.307 

1.0393 

8.537 

8500 

36.115(3) 

0.9643(1) 

0.0357(1) 

0.4796(5) 

4.682(2) 

12.1932 

14.352 

1.0863 

8.536 

8550 

36.119(3) 

0.9632(1) 

0.0368(1) 

0.4800(8) 

4.687(2) 

12.3059 

14.455 

1.1950 

8.539 
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Table  3.  Constant-pressure  RxMC  simulations  of  shocked  liquid  N2  (continued). 


T 

<P> 

Mole  Fraction3 

<V> 

<lfoaf> 

H°b 

Ec 

Hg“ 

D e 

(K) 

(GPa) 

<x(N2)> 

<  x(N)> 

(cm3/g) 

(kJ/g) 

(kJ/g) 

(kJ/g) 

(kJ/g) 

(km/s) 

Pimp  =  47.0  GPa 

9350 

47.152(5) 

0.9260(1) 

0.0741(1) 

0.4351(8) 

5.809(3) 

14.9046 

17.939 

-0.5608 

9.481 

9450 

47.182(2) 

0.9221(2) 

0.0779(2) 

0.4352(6) 

5.812(2) 

15.2050 

18.212 

-0.2958 

9.485 

9550 

47.148(2) 

0.9183(2) 

0.0816(2) 

0.4355(3) 

5.805(2) 

15.5024 

18.473 

-0.0142 

9.484 

9575 

47.130(2) 

0.9179(1) 

0.0821(1) 

0.4355(3) 

5.799(2) 

15.5580 

18.515 

0.0353 

9.482 

9600 

47.134(3) 

0.9166(2) 

0.0834(2) 

0.4357(6) 

5.804(2) 

15.6475 

18.602 

0.1249 

9.483 

Pimp  =  52.6  GPa 

9500 

52.740(4) 

0.9101(2) 

0.0899(2) 

0.4160(6) 

6.336(3) 

15.7319 

19.248 

-1.9974 

9.911 

9750 

52.737(3) 

0.9002(2) 

0.0998(2) 

0.4162(4) 

6.325(2) 

16.5027 

19.934 

-1.3053 

9.912 

10000 

52.754(3) 

0.8888(2) 

0.1112(2) 

0.4162(4) 

6.314(3) 

17.3362 

20.683 

-0.5620 

9.914 

10250 

52.755(6) 

0.8774(2) 

0.1226(2) 

0.4164(8) 

6.304(4) 

18.1791 

21.441 

0.2006 

9.915 

Pimp  =  60.4  GPa 

10400 

60.527(3) 

0.8498(3) 

0.1502(3) 

0.3925(5) 

6.989(2) 

19.4771 

23.379 

-1.7771 

10.469 

10500 

60.521(4) 

0.8441(2) 

0.1559(2) 

0.3924(5) 

6.979(2) 

19.8615 

23.724 

-1.4344 

10.468 

10600 

60.563(7) 

0.8388(3) 

0.1612(3) 

0.3920(8) 

6.969(4) 

20.2337 

24.057 

-1.1294 

10.469 

10900 

60.551(4) 

0.8231(2) 

0.1769(2) 

0.3916(7) 

6.939(2) 

21.3514 

25.056 

-0.1370 

10.466 

10950 

60.506(3) 

0.8201(3) 

0.1799(3) 

0.3916(6) 

6.928(4) 

21.5531 

25.231 

0.0583 

10.462 

11000 

60.544(4) 

0.8170(4) 

0.1830(4) 

0.3915(8) 

6.929(2) 

21.7593 

25.424 

0.2306 

10.465 

Pimp  =  81.1  GPa 

12400 

81.093(5) 

0.6487(4) 

0.3513(4) 

0.3352(5) 

8.275(4) 

30.8436 

35.438 

-0.7344 

11.728 

12500 

81.075(7) 

0.6416(5) 

0.3584(5) 

0.3347(7) 

8.253(3) 

31.3115 

35.854 

-0.3281 

11.724 

12600 

81.120(8) 

0.6351(4) 

0.3649(4) 

0.3343(8) 

8.240(4) 

31.7590 

36.259 

0.0387 

11.724 

Pimp  =  91.5  GPa 

13000 

91.191(9) 

0.5584(6) 

0.4416(6) 

0.3109(7) 

8.775(5) 

35.5438 

40.461 

-1.3754 

12.274 

13250 

91.118(9) 

0.5424(5) 

0.4576(5) 

0.3101(8) 

8.725(4) 

36.6715 

41.464 

-0.3756 

12.263 

13500 

91.092(11) 

0.5267(7) 

0.4733(7) 

0.3092(8) 

8.671(6) 

37.7904 

42.454 

0.5857 

12.255 

15500 

90.482(12) 

0.4170(6) 

0.5830(6) 

0.3045(11) 

8.242(3) 

46.2912 

49.933 

8.1337 

12.183 

16000 

90.418(14) 

0.3917(6) 

0.6083(6) 

0.3036(12) 

8.147(4) 

48.3805 

51.779 

9.9691 

12.173 

16500 

90.064(18) 

0.3599(8) 

0.6401(8) 

0.3022(17) 

7.991(4) 

50.7644 

53.858 

12.1530 

12.141 

Notes:  ‘<  >’  indicates  ensemble  averages  detennined  from  the  simulation. 

Uncertainty  in  units  of  the  last  decimal  digit  is  given  in  parentheses:  3.010(2)  means  3.010  ±  0.002. 

“Mole  fraction  of  N2,  so  x(N2)  =NN2/Ntotal  and  x(N)  =  lA  NN/Ntotal,  where  N  total  =  3375. 

bFrom  equation  11. 

cFrom  equation  14. 

dFrom  equation  8. 

eFrom  equation  7. 

experimental  data  is  given  in  table  4.  In  table  4,  the  uncertainties  for  the  shock  Hugoniot 
properties  measured  by  experiment  are  given  in  parentheses,  while  the  uncertainties  in  the  RxMC 
calculations  can  be  estimated  from  the  R-square  value  of  the  functional  fit  of  data  given  in 
table  3.  Typical  R-square  values  for  the  predicted  temperatures  and  specific  volumes  are 
0.97-0.99.  Plots  of  the  shock  Hugoniot  pressure  vs.  the  specific  volume  and  the  shock  wave 
velocity  are  given  in  figures  1  and  2,  respectively,  for  both  the  reactive  and  nonreactive  models. 
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Table  4.  Shock  Hugoniot  states  of  liquid  nitrogen.3  RxMC  results  are  for  the  reactive  model  discussed  in  the 
text.  Experimental  data  is  taken  from  Nellis  et  al.  (55),  except  for  those  noted. 


p 

(GPa) 

V 

(cm3/mol  N2) 

T 

(K) 

D 

(km/s) 

Exp. 

RxMCb 

Exp. 

RxMC 

Exp. 

RxMC 

Exp. 

RxMC 

2.96° 

2.96 

21.7 

21.36 

— 

536.2 

3.14 

3.087 

4.74° 

4.74 

20.1 

19.82 

— 

883.9 

3.74 

3.703 

10.1c 

10.0 

17.3 

17.41 

— 

2008.4 

5.00 

4.993 

18.1(0.5) 

18.1 

15.34(0.5) 

15.57 

4300(200) 

3912.4 

6.34 

6.384 

29.9(0.5) 

29.9 

14.26(0.5) 

14.05 

7300(250) 

6778.1 

7.93 

7.895 

36.0(0.4) 

36.0 

13.41(0.4) 

13.35 

8750(300) 

7963.0 

8.52 

8.519 

47.0(0.5) 

47.0 

11.83(0.5) 

12.20 

8900(600) 

9557.7 

9.40 

9.483 

52.6(0.5) 

52.6 

11.13(0.5) 

11.66 

11100(800) 

10185.4 

9.79 

9.914 

60.4(0.7) 

60.4 

10.31(0.7) 

10.97 

12000(850) 

10935.2 

10.31 

10.465 

81.1(1.5) 

81.1 

9.40(1.5) 

9.366 

14500(1000) 

12588.9 

11.73 

11.724 

“Uncertainties  in  experimental  data  (where  available)  are  given  in  parentheses. 
bPressure  imposed  in  the  constant-pressure  version  of  the  RxMC  method. 
“Taken  from  Zubarev  and  Telegin  (54). 


volume  [cm3/mol  of  N2] 


Figure  1.  Shock  Hugoniot  of  liquid  N2.  Calculated  values  from  RxMC  simulations  using  a 
reactive  (o)  and  nonreactive  (□)  model  are  compared  with  experimental  data  (A) 
(54,  55).  The  shock  pressure  is  plotted  vs.  the  molar  volume  of  N2. 
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Figure  2.  Shock  Hugoniot  of  liquid  N2.  Calculated  values  from  RxMC  simulations  (reactive 
model:  [o];  nonreactive  model  [□])  are  compared  with  experimental  data  (A)  (54, 

55).  The  shock  pressure  is  plotted  vs.  the  shock  wave  velocity. 

We  found  excellent  agreement  between  the  RxMC  calculations  using  the  reactive  model  and  the 
experimental  measurements  for  most  of  the  pressures  considered.  As  the  pressure  is  increased 
along  the  Hugoniot,  the  system  becomes  a  partially  dissociated  fluid  containing  a  mixture  of  N2 
molecules  and  N  atoms.  This  behavior  is  reflected  in  the  species  mole  fractions  plot  of  figure  3. 
Values  of  the  mole  fractions  shown  in  figure  3  were  determined  by  interpolating  the  data 
reported  in  table  3  to  THg  using  a  quadratic  function.  It  is  evident  from  figure  1  that  the 
nonreactive  model  fails  at  high  pressures  where  the  dissociation  is  not  negligible  (11). 

The  experimental  data  appears  to  exhibit  a  softening  of  the  Hugoniot  curve  near  55  GPa  while 
the  RxMC  simulations  do  not  predict  this  behavior.  Experimental  errors  are  increasing  in  this 
region;  therefore,  whether  the  discrepancy  is  due  to  experimental  uncertainty  or  an  inaccurate 
model  cannot  be  conclusively  established  (11).  Interestingly,  however,  recent  density-functional 
theory  (DFT)  calculations  (31)  predict  the  softening  behavior.  Pair  correlation  calculations  in 
the  work  of  Kress  et  al.  (31)  indicate  that  the  system  contains  a  small  fraction  of  clusters  larger 
than  dimers  in  the  partially  dissociated  region;  however,  these  larger  clusters  are  not  long-lived. 
The  dissociated  system  may  contain  N„  molecules  bound  by  single  (N-N)  or  single  and  double 
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P  [GPa] 


Figure  3.  Species  mole  fractions  (N2:  [♦];  N:  [A])  along  the  Flugoniot  curve  determined  from 
RxMC  simulations  of  the  N2  dissociation  reaction. 

N  =  N)  bonds.  Although  refining  the  N2  model  is  outside  the  scope  of  the  present  work,  such 
products  and  their  accompanying  reactions  could  be  included  into  the  RxMC  simulation 
scenario,  e.g.,  N  =  NON  =  N  or  N  =  NON-N,  which  may  make  RxMC  calculations  more 
compatible  with  the  experimental  measurements. 

As  a  matter  of  curiosity,  a  point  along  the  Hugoniot  curve  at  P  =  91.5  GPa  was  calculated. 
Although  experimental  data  are  not  presently  available  at  this  pressure,  the  recent  DFT 
calculations  of  Kress  and  coworkers  (31)  predict  considerably  different  behavior  in  this  pressure 
regime.  The  model  used  in  the  present  work  predicts  softening  behavior  (v  =  8.68  cm  /mole  N2), 
as  does  the  work  of  Fried  and  Howard  (11),  while  the  DFT  calculations  appear  to  be  approaching 
a  maximum  compression  in  this  region  (v  =  9.34  cm3/mole  N2  atP  =  91.5  GPa  [3f]).  According 
to  our  calculations,  such  a  system  would  contain  a  nearly  equimolar  mixture  of  N2  molecules  and 
dissociated  N  atoms  (see  figure  3). 
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5.  Discussion 


We  have  demonstrated  the  effectiveness  of  using  the  RxMC  simulation  method  for  determining 
the  shock  properties  of  materials.  We  found  the  RxMC  calculations  to  be  in  excellent  agreement 
with  the  available  experimental  data  for  the  N2  dissociation  reaction.  These  results  illustrate  the 
utility  of  the  method  for  predicting  the  shock  Hugoniot  for  mixtures  in  which  species 
concentrations  are  not  known  and  in  the  absence  of  interaction  potentials  that  simulate  bond 
breakage  and  formation.  In  the  next  report  of  this  series,  we  consider  an  application  of  the 
RxMC  method  to  the  nitric  oxide  decomposition  reaction  (2NOON2+O2).  Several  possible 
extensions  of  the  RxMC  are  also  discussed  in  the  subsequent  report  {34). 
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B  HOMAN 
A  HORST 
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S  HOWARD 
PKASTE 
G  KATULKA 
T  KOGLER 
A  KOTLAR 
C  LEVERITT 
RLIEB 
D  LYON 
K  MCNESBY 
M  MCQUAID 
M  MILLER 
A  MIZIOLEK 
J  MORRIS 
J  NEWBERRY 
J  NEWILL 
M  NUSCA  (6  CPS) 

W  OBERLE 
R  PESCE-RODRIGUEZ 
P  PLOSTINS 
S  PIRIANO 
G  REEVES 
B  RICE 
JSAHU 
R  SAUSA 
S  SILTON 
E  SCHMIDT 
J  SCHMIDT 
P  WEINACHT 
D  WILKERSON 
A  WILLIAMS 
M  ZOLTOSKI 
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1  ACADEMY  OF  SCIENCES 
OF  THE  CZECH  REPUBLIC 
E  HALA  LABORATORY 
OF  THERMODYNAMICS 
M  LISAL 

165  02  PRAGUE  6-SUCHDOL 
CZECH  REPUBLIC 
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